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Abstract 

We introduce, for the first time, a new class of Birnbaum-Saunders nonlinear 
regression models potentially useful in lifetime data analysis. The class generalizes 
the regression model described by Rieck and Nedelman [1991, A log-linear model 
for the Birnbaum-Saunders distribution, Technometrics, 33, 51-60]. We discuss 
maximum likelihood estimation for the parameters of the model, and derive closed- 
form expressions for the second-order biases of these estimates. Our formulae are 
easily computed as ordinary linear regressions and are then used to define bias 
corrected maximum likelihood estimates. Some simulation results show that the 
bias correction scheme yields nearly unbiased estimates without increasing the mean 
squared errors. We also give an application to a real fatigue data set. 

Key words: Bias correction, Birnbaum-Saunders distribution, maximum 
likelihood estimation, nonlinear regression. 



1 Introduction 



Different regression models have been proposed for lifetime data such as those 
based on the gamma, lognormal and Weibull distributions. These models typi- 
cally provide a satisfactory fit in the middle portion of the data, but very often 
fail to deliver a good fit at the tails, where only a few observations are gene- 
rally available. The family of distributions proposed by Birnbaum and Saun- 
ders (1969) can also be used to model lifetime data and it is widely applicable 
to model failure times of fatiguing materials. This family has the appealing 
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feature of providing satisfactory tail fitting. This family of distributions was 
originally obtained from a model for which failure follows from the develop- 
ment and growth of a dominant crack. It was later derived by Desmond (1985) 
using a biological model which followed from relaxing some of the assumptions 
originally made by Birnbaum and Saunders (1969). 

The random variable T is said to be Birnbaum-Saunders distributed with 
parameters a, 77 > 0, say B-S{a,ri), if its cumulative distribution function 
(cdf) is given by 



where $(•) is the standard normal distribution function and a and rj are shape 

and scale parameters, respectively. It is easy to show that rj is the median of 
the distribution: Fr(ry) = ^{0) = 1/2. For any k > 0, then kT ~ B-S{a, krj). 

McCarter (1999) considered parameter estimation under type II data censo- 
ring for the B-S{a, rj) distribution. Lemonte et al. (2007) derived the second- 
order biases of the maximum likelihood estimates (MLEs) of a and 1], and 
obtained a corrected likelihood ratio statistic for testing the parameter a. 
Lemonte et al. (2008) proposed several bootstrap bias corrected estimates of 
a and 77. Further details on the Birnbaum-Saunders distribution can be found 
in Johnson et al. (1995). 

Rieck and Nedelman (1991) proposed a log-linear regression model based on 
the Birnbaum-Saunders distribution. They showed that if T ~ B-S{a,7]), 
then Y — log(T) is sinh- normal distributed, say Y ~ SM{a, ji, a), with shape, 
location and scale parameters given by a, = log(ry) and cr = 2, respectively. 
Their model has been widely used as an alternative model to the gamma, 
lognormal and WeibuU regression models; see Rieck and Nedelman (1991, § 7). 
Diagnostic tools for the Birnbaum-Saunders regression model were developed 
by Galea et al. (2004), Leiva et al. (2007) and Xie and Wei (2007), and the 
Bayesian inference was considered by Tisionas (2001). 

In this paper we propose a class of Birnbaum-Saunders nonlinear regression 
models which generalizes the regression model introduced by Rieck and Nedel- 
man (1991). We discuss maximum likelihood estimation of the regression pa- 
rameters and obtain the Fisher information matrix. As is well known, however, 
the MLEs, although consistent, are typically biased in finite samples. In order 
to overcome this shortcoming, we derive a closed-form expression for the bias 
of the MLE in these models which is used to define a bias corrected estimate. 

Bias adjustment has been extensively studied in the statistical literature. In 
fact. Cook et al. (1986) proposed bias correction in normal nonlinear models. 
Young and Bakir (1987) obtained bias corrected estimates for a generahzed 
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log- gamma regression model. Cordeiro and McCuUagh (1991) gave general ma- 
trix formulae for bias correction in generalized linear models, whereas Paula 
(1992) derive the second-order biases in exponential family nonlinear models. 
Cordeiro et al. (2000) obtained bias correction for symmetric nonhnear re- 
gression models. More recently, Vasconcellos and Cribari-Neto (2005) calcu- 
late the biases of the MLEs in a new class of beta regression. Cordeiro and 
Demetrio (2008) propose formulae for the second-order biases of the maxi- 
mum quasi-likelihood estimates, whereas Cordeiro and Toyama (2008) derive 
the second-order biases in generalized nonhnear models with dispersion co- 
variates. 

The rest of the paper is as follows. Section 2 introduces the class of Birnbaum- 
Saunders nonlinear regression models and discusses maximum likelihood esti- 
mation. Using general results from Cox and Snell (1968), we derive in Section 3 
the second-order biases of the MLEs of the nonlinear parameters in our class 
of models and define bias corrected estimates. Some special models are consi- 
dered in Section 4. Simulation results are presented and discussed in Section 5 
for two nonlinear regression models. We show that the bias corrected estimates 
are nearly unbiased with mean squared errors very close to the corresponding 
ones of the uncorrected estimates. Section 6 gives an application of the pro- 
posed regression model to a real fatigue data set, which provides a better fit 
at the tail of the data. Finally, Section 7 concludes the paper. 



2 Model specification 



Let T ~ B-S{a,r]). The density function of F = log(T) ~ SJ\f{a, iJ,,a) has 
the form (Rieck and Nedelman, 1991) 

7r(y; a, cr) = ^"^^ (^~^) ^^^{~^^^^^^ (^V^) }' ^ ^ 

This distribution has a number of interesting properties (Rieck, 1989): (i) It 
is symmetric around the location parameter (ii) It is unimodal for a < 2 
and bimodal for a > 2; (iii) The mean and variance of Y are K{Y) = fi and 
Var(F) = a'^w{a), respectively. There is no closed-form expression for w{a), 
but Rieck (1989) obtained asymptotic approximations for both small and large 
values of a; (iv) If Ya ~ SJ\f{a, fj,, a), then Sa — 2{Ya — A*)/ (aa) converges in 
distribution to the standard normal distribution when a — > 0. 

We define the nonlinear regression model 

Vi ^ Mxi, f3) + Ei, i = l,...,n, (1) 

where is the logarithm of the ith observed lifetime, ajj is an m x 1 vector 
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of known explanatory variables associated with the ith observable response 
Hi, /3 = . . . ^I3.p)^ is a vector of unknown nonlinear parameters, and Ei ~ 
SJ\f{a, 0, 2). We assume a nonlinear structure for the location parameter /Xj in 
model (1), say /Xj = fi{Xi; /3), where fi is assumed to be a known and twice 
continuously differentiable function with respect to (3. For the linear regression 
Hi — xj 13, the model (1) reduces to Rieck and Nedelman's (1991) model. 

The log-likelihood function for the vector parameter d = (/3^,q;)^ from a 
random sample y — (yi, . . . , y^)^ obtained from (1), except for constants, can 
be expressed as 

n 1 

m^j:^ogiCa)-::T.Cl (2) 

1=1 ^ i=i 

where = ^ii{e) = 2a-^cosh{[yi-fii]/2),^i2 = ^,2(0) = 2a"i sinh([|/i-/ii]/2) 
for i = 1, . . . , n. The function £{6) is assumed to be regular (Cox and Hinkley, 
1974, Ch. 9) with respect to all (3 and a derivatives up to third order. Further, 
the n X p local matrix D — D(f3) — dn/dfS of partial derivatives of /x = 
. . . , Hn)^ with respect to (3 is assumed to be of full rank, i.e., rank(Z)) = p 
for all (3. The nonlinear predictors embedded in an infinite 

sequence of m x 1 vectors that must satisfy these regularity conditions for 
the asymptotics to be vahd. Under these assumptions, the MLEs have good 
asymptotic properties such as consistency, sufficiency and normality. 

The derivatives with respect to the components of (3 and a are denoted by: 
Ur = de{e)/d(3r, Ua = di{d)/da, U^s = dH{d)/d(3rd(3s: Ura = dH{0)/d(3rda, 
Urst = d^e{e)/dPrdPsdPt, Ursa = dH{e) /dprOPsda, etc. Further, we use the 
following notation for joint cumulants of log-likelihood derivatives: Krs = 

K{Urs), K,r,a = ^{UrUa), K-rst = '^(Urst), Ctc. Let = dKrs/d/3t, CtC. All 

«;'s and their derivatives are assumed to be of order 0{n). Also, we adopt the 
notation dir — diii/df3r and girs = d'^/ii/dPrdPs for the first and second partial 
derivatives of /li with respect to the elements of /3. 

It is easy to see by differentiating (2) that 

Ur — -^y2 ( ^ilCi2 ~ 7^ V Ua — \ ^^i2: 

^ i=l \ ^ilj ^ ^ i=l 

1 " n 3 " 

Ura — y,dir^ii$,i2 and Uaa — 5- / ,^i2- 

The score function for (3 is t/g — |-D^s, where s — s{d) is an n-vector whose 
ith element is equal to ^jiCi2 ~ Ci2/ Cu- 
lt is well-known that, under general regularity conditions (Cox and Hin- 
kley, 1974, Ch. 9), the MLEs are consistent, asymptotically efficient and 
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asymptotically normal. Let 6 = [p' , S)^ be the MLE oi — (/3^, a)^. We 
can write ~ A/^+i(0, li'g ^) for n large, where ~ denotes approximately 
distributed, Kg is the block-diagonal Fisher information matrix given by 
Ke = dia,g{K K,a,a} , is its inverse, = tjji{a){D^ D)/ A is the in- 

formation matrix for /3 and K.a,a = '^njo? is the information for a. Also, 

where erf (•) is the error function given by 

erf (x) = ^ / e"* dt. 

\Jt^ Jo 

Details on erf (■) can be found in Gradshteyn and Ryzhik (2007). Since Kq 
is block-diagonal, the vector (3 and the scalar a are globally orthogonal (Cox 
and Reid, 1987) and (3 and a are asymptotically independent. It can be shown 
(Rieck, 1989) that il)i{a) k, 1 -f for a small and ■?/'i(a) pa 2 for a large. 

The MLE d satisfies p+1 equations Ur = [/« = for the components of f3 and 
a. The Fisher scoring method can be used to estimate /3 and a simultaneously 
by iteratively solving the equations 

p{m+l) _ p{m) _|_ ^]j{m)T ■Q{m)-j-l j~f{rri)T ^{m) ^ 

where C^""^ = 2s('"VV'i(a^'"^) and C?"^ = Ei=i ilt^ /n for m = 0, 1, 2, . . .. 

The above equations show that any software with a weighted linear regression 
routine can be used to calculate the MLEs of (3 and a iteratively. Initial 
approximations (3^^^ and a*^"-* for the iterative algorithm are used to evaluate 
£)(o) ^ ^(0) and ^2°^ from which these equations can be used to obtain the next 
estimates (3'^^^ and a^^\ These new values can update £>, C, and ^2 and so the 
iterations continue until convergence is achieved. 



3 Biases of estimates of j3 and a 



We now obtain some joint cumulants of log-likelihood derivatives and their 
derivatives: 

__V;i(a)" _ _ __2n _10n 

4: ~[ a"^ 
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Ma) ^^ , _ (2 + Q^) ^ ^ ^ 

^rst — , / XQirs^it ~r Qirt^is ~l" 9ist^ir) i ^rsa q /.^ir^isi 

(t) _ ^i(Q') y-/ J. , J. N (a) _ (s) _ Q ^ (a) _ ^ 

^ i=l ^ 

Let B{j3a) and B{a) be the biases of (a = 1, . . . , p) and 3, respectively. 
The use of Cox and Snell's (1968) formula to obtain these biases is greatly sim- 
plified, since (3 and a are globally orthogonal and the cumulants corresponding 
to the parameters in {3 are invariant under permutation of these parameters. 
Prom now on we use Einstein summation convention with the indices varying 
over the corresponding parameters. We have 



(3) 



and 

B{a) = {K'^n' - l^aaa^ + Z^"'" E' ^''^ (^^i^ ' ^^atu^ , (4) 

where m:''* is the (r, s)th element of the inverse of the information matrix 
for f3, k"'" = K~]^ and Yl denotes the summation over all combinations of 
parameters (3i, (32, ■ ■ ■ , Pp. 

First, we consider equation (3) from which we readily have that the second 
sum is zero since Ugaa — i^^^ = 0. It follows that 

B0a) = -^^^ E' Z^''''/^*'" Y.{gisudu - Qistdiu + Qitudis). 
s,t,u i=l 

By rearranging the summation terms we obtain 

B0a) = -^E E'^^'^'^^^^E''^*'"^^*"- 

° 1=1 s t,u 

Let dj {1 X p) and gj (1 x p^) be vectors containing the first and second 
partial derivatives of the mean /ij with respect to the /9's. We can write the 
above equation in matrix notation as 



° i=i 

where is the ath row of the pxp identity matrix and vec(-) is the operator 
which transforms a matrix into a vector by stacking the columns of the matrix 
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one underneath the other. It is straightforward to check that 

where D = dfi/d(3 = (di, . . . , dn)^ and G = d^n/d(3^dl3 = {g^ . . .^gnf 
are n x p and n x matrices of the first and second partial derivatives of the 
mean vector /x with respect to /3, respectively. The bias vector B{f3) of 
^ can then be written as 

B0) = (D^D) 'D^d, (5) 
where d is an n x 1 vector defined as cZ = — {2/V'i(Q;)}Gvec{(£>^£>)~^}. 
We now calculate the bias of a. Using (4), we obtain 

« (2 + «') ^ t«j ^ « (2 + «') ^^^T,v-l^ 



4n 4Q;n 



-tr(Z)K^^D'), 



where tr(-) denotes the trace operator. Now, making use of the fact that 
ix^DK'p^ ) = 4:p/ipi{a), we can rewrite the bias of a as 

1 ( f 2 + a'^ \ a] 
B{a) = —{p[^— -]+-}. 6 



Equations (5) and (6) represent the main results of the paper. The bias vector 
B{f5) can be obtained from a simple ordinary least-squares regression of d 
on the columns of D. It depends on the nonlinearity of the regression func- 
tion / and the parameter a. The bias vector B((3) will be small when d is 
orthogonal to the columns of D. Also, it can be large when ^l)l{a) and n are 
both small. Equation (5) is easily handled algebraically for any type of non- 
linear regression, since it involves simple operations on matrices and vectors. 
For special models with closed-form information matrix for /3, it is possible 
to obtain closed-form expressions for B0). For linear models, the matrix G 
and the vector d vanish and hence B0) — 0, which is in agreement with the 
result due to Rieck and Nedelman (1991, p. 54) that the MLEs arc unbiased 
to order n~^. Expression (6) depends directly on the nonlinear structure of 
the regression model only through the rank p of D. It shows that the bias is 
always a linear function of the dimension p 13. 

In the right-hand sides of expressions (5) and (6), which are both of order , 
consistent estimates of the parameters (3 and a can be inserted to define bias 
corrected estimates f3 = j3 — B{f3) and a = a — B{a), where B{f3) and B{a) 
are the values of B((3) and -B(a), respectively, at 6 — The bias 
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corrected estimates f3 and a are expected to have better sampling properties 
than the classical MLEs f3 and a. In fact, we present some simulations in 
Section 5 to show that (3 and a have smaller biases than their corresponding 
uncorrected estimates, thus suggesting that these bias corrections have the 
effect of shrinking the adjusted estimates toward to the true parameter values. 
However, we can not say that the bias corrected estimates offer always some 
improvement over the MLEs, since they can have mean squared errors larger. 

It is worth emphasizing that there are other methods to obtain bias corrected 
estimates. In regular parametric problems. Firth (1993) developed the so- 
called "preventive" method, which also allows for the removal of the second- 
order bias. His method consists of modifying the original score function to 
remove the first-order term from the asymptotic bias of these estimates. In 
exponential families with canonical parameterizations, his correction scheme 
consists in penalizing the likelihood by the Jeffreys invariant priors. This is a 
preventive approach to bias adjustment which has its merits, but the connec- 
tions between our results and his work are not pursued in this paper since they 
could be developed in future research. Additionally, it should be mentioned 
that it is possible to avoid cumbersome and tedious algebra on cumulant cal- 
culations by using Efron's bootstrap (Efron and Tibshirani, 1993). We use the 
analytical approach here since this leads to a nice formula. Moreover, the ap- 
plication of the analytical bias approximation seems to generally be the most 
feasible procedure to use and it continues to receive attention in the literature. 

We now calculate the second-order bias B{jii) of the MLE fxi of the ith mean 
— fi{xi] 13). We can easily show by Taylor series expansion that 

S(/2,) = d]B{^) + itr{M,Cov(^)}, 

where Mi is a p x p matrix of second partial derivatives d^/ii/dPrdPs (for 
r, s = l,...,p), Cov(;9) = K^^ is the asymptotic covariance matrix of ^ 
and the vectors di and B{j3) were mentioned previously. All quantities in the 
above equation should be evaluated at /3. 

The asymptotic variance of jii can also be expressed explicitly in terms of the 
covariance of /9 by 

Var(/20 = tT{{didJ )Coy0)}. 



4 Special models 

Equation (5) is easily handled algebraically for any type of nonlinear model, 
since it involves simple operations on matrices and vectors. This equation, 
in conjunction with a computer algebra system such as MAPLE (Abell and 
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Braselton, 1994) will compute B{f3) algebraically with minimal effort. In par- 
ticular, (5) may simplify considerably if the number of nonlinear parameters 
is small. Moreover, for any nonlinear special model, we can calculate the bias 
B(f3) numerically via a software with numerical linear algebra facilities such 
as Ox (Doornik, 2001) and R (R Development Core Team, 2008). 

First, we consider a nonlinear regression model which depends on a single 
nonlinear parameter. Equation (5) gives 

V'i(a) «i 

where Ki = E2=i(dfi/d(3f and K2 = E7=i{df^/d[3){d'^ fi/d(3^). The constants 
Ki and K2 are evaluated at f3 and a to yield B{f3) and the corrected estimate 
j3 = (3 — B{j3). For example, the simple exponential model /j = exp(/3a;j) leads 
to = Er=ia;- exp(2/3a;i) and K2 = EjLi exp(2/3a;i). 

As a second example, we consider a partially nonlinear regression model de- 
fined by 

tx^ZX + rig{j), (7) 

where .Z is a known n x (p — 2) matrix of full rank, g{'y) is an n x 1 vector, 
f3 = (A""", 7/, 7)^, A = (Ai, . . . , Ap_2)''' and t] and 7 are scalar parameters. 
This class of models occurs very often in statistical modeling; see Cook et 
al. (1986) and Cordeiro et al. (2000). For example, = Ai^ii -|- A2^2 + '7exp(7a;) 
(Gallant, 1975), n ~ X — 77log(a;i -|- 70:2) (Darby and Ellis, 1976) and fi — 
A-|-77log(xi/(7-|-X2)) (Stone, 1980). Ratkowsky (1983, Ch. 5) discusses several 
models of the form (7) which include the asymptotic regression and WeibuU- 
type models given hy /i — X — rjj^ and /i — X — 'r]exp{—jx), respectively. 

The n X p local model matrix D takes the form D — [Z,g{'y),r]{dg{'y)/d'y)] 
and, after some algebra, we can obtain from (5) 



B{f3) 



-Coy{fj, j)Tp + ^Var(7)5p 



(8) 



where Tp is a p x 1 vector with a one in the last position and zeros else- 
where, dp = {D^D)^^D^{d'^g{'y)/d'y'^) is simply the set of coefficients from 
the ordinary regression of the vector d^gf(7)/d7^ on the matrix £>, and Var(7) 
and Cov(f7, 7) are the large-sample second moments obtained from the ap- 
propriate elements of the asymptotic covariance matrix Cov(;9) = K^-^ = 
{4/ilJi{a)){D^ D)^^ . It is clear from (8) that B{(3) does not depend explicitly 
on the linear parameters in A and it is proportional to A/ipi{a). Further, the 
covariance term Cov(f/, 7) contributes only to the bias of 7. 
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5 Numerical results 



We now use Monte Carlo simulation to evaluate the finite-sample performance 
of the MLEs of the parameters and of their corrected versions in two nonlinear 
regression models. The MLEs of the parameters were obtained by maximizing 
the log-hkehhood function using the BFGS quasi-Newton method with ana- 
lytical derivatives. This method is generally regarded as the best-performing 
nonlinear optimization method (Mittelhammer et al., 2000, p. 199). The co- 
variate values were selected as random draws from the uniform U{0, 1) distri- 
bution and for fixed n those values were kept constant throughout the exper- 
iment. Also, the number of Monte Carlo replications was 10,000. All simula- 
tions were performed using the Ox matrix programming language (Doornik, 
2001).^ 

In order to analyze the performance of the estimates, we computed, for each 
sample size and for each estimate: the relative bias (the relative bias of an 
estimate 9, defined as {E{0) — 9}/9, is obtained by estimating E{9) by Monte 
Carlo) and the root mean square error (VMSE), where MSE is the estimated 
mean square error from the 10,000 Monte Carlo replications. 

First, we consider the nonlinear regression model 

A«j = >^iZii + \2Zi2 77exp(7Xj), 

where Si ~ SJ\f{a, 0, 2) for i = 1, ... , n. The sample sizes were n = 15, 30 and 
45. Without loss of generality, the true values of the regression parameters 
were taken as Ai = 4, A2 = 5, 77 = 3, 7 = 1.5 and a — 0.5 and 1.5. 

Table 1 gives the relative biases of both uncorrected and corrected estimates to 
show that the bias corrected estimates are much closer to the true parameters 
than the unadjusted estimates. For instance, when n — 15 and a — 1.5, 
the average of the estimated relative biases for the estimates of the model 
parameters is —0.03244, whereas the average of the estimated relative biases 
for the corrected estimates is —0.0083. Hence, the average bias (in absolute 
value) of the MLEs is almost four times greater than the average bias of the 
corrected estimates. This fact suggests that the second-order bias of the MLEs 
should not be ignored in samples of small to moderate size since they can be 
non-negligible. The figures in Table 2 show that the root mean squared errors 
of the uncorrected and corrected estimates are very close. Hence, the figures 
in both tables suggest that the corrected estimates have good properties. 

When the parameter a increases, the finite-sample performance of the MLEs 



^ Ox is freely distributed for academic purposes and available at 
http://www.doornik.com. 
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Table 1 

Relative biases of the uncorrected and corrected estimates. 

an Ai A2 ?7 7 a 

0.5 15 MLE 0.0006 -0.0013 0.0011 0.0020 -0.1691 

BCE 0.0007 -0.0011 0.0001 0.0008 -0.0395 

30 MLE 0.0001 -0.0013 0.0013 0.0009 -0.0811 

BCE 0.0002 -0.0012 0.0007 -0.0001 -0.0092 

45 MLE 0.0003 -0.0012 0.0007 0.0008 -0.0537 

BCE 0.0003 -0.0011 0.0003 0.0001 -0.0042 

1.5 15 MLE -0.0068 -0.0083 0.0248 0.0197 -0.1916 

BCE -0.0055 -0.0046 0.0113 0.0056 -0.0481 

30 MLE -0.0016 -0.0034 0.0079 0.0078 -0.0933 

BCE -0.0011 -0.0018 0.0027 0.0012 -0.0116 

45 MLE -0.0028 -0.0027 0.0052 0.0026 -0.0614 

BCE -0.0023 -0.0018 0.0023 -0.0005 -0.0048 

BCE: bias corrected estimate. 

Table 2 

Root mean squared errors of the uncorrected and corrected est imates . 
an Ai A2 ?? 7 a 

0.5 15 MLE 0.4093 0.4920 0.2707 0.0924 0.1234 

BCE 0.4093 0.4921 0.2709 0.0922 0.1067 

30 MLE 0.3006 0.3806 0.2113 0.0688 0.0763 

BCE 0.3006 0.3806 0.2114 0.0686 0.0702 

45 MLE 0.2434 0.2874 0.1768 0.0567 0.0590 

BCE 0.2434 0.2874 0.1769 0.0566 0.0555 

1.5 15 MLE 1.6302 1.1230 0.9756 0.3235 0.3938 

BCE 1.6333 1.1274 0.9819 0.3152 0.3315 

30 MLE 0.9684 0.7003 0.5785 0.1931 0.2399 

BCE 0.9693 0.7011 0.5807 0.1908 0.2155 

45 MLE 0.6505 0.5575 0.3895 0.1318 0.1837 

BCE 0.6507 0.5577 0.3901 0.1311 0.1700 

BCE: bias corrected estimate. 
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deteriorates (see Tables 1 and 2). For instance, when n — 15, the relative biases 
of 7 (MLE) and 7 (BCE) were 0.0020 and 0.0008 (for a = 0.5) and 0.0197 
and 0.0056 (for a = 1.5), which indicate an increase in the relative biases of 
nearly 10 and 7 times, respectively. Also, the root mean squared errors in the 
same order were 0.0924 and 0.0922 (for a = 0.5) and 0.3235 and 0.3152 (for 
1.5). 

Next, we consider the very known Michaelis-Menton model, which is very use- 
ful for estimating growth curves, where it is common for the response to ap- 
proach an asymptote as the stimulus increases. The Michaelis-Menton model 
(Ratkowsky, 1983) provides an hyperbolic form for /Xj against Xi given by 

Hi = — ■ , z = l,2, 

where the curve has an asymptote at = r^. Here, the sample sizes were 
n — 20, 30, 40 and 50. Also, the true values of the regression parameters were 
taken as 77 = 3 and 7 = 0.5, with a — 0.5. 

Table 3 gives the relative biases and root mean squared errors of the uncor- 
rected and corrected estimates. The figures in this table reveal that the MLEs 
of the parameters can be substantially biased, even when n — 50, and that 
the bias correction is very effective. In terms of MSE, the adjusted estimates 
are slightly better than the ordinary MLEs. 

Table 3 

Relative biases and root mean squared errors of the uncorrected and corrected 
estimates; a = 0.5 and different sample sizes. 



Relative Bias 


VMSE 


n rj J a 


7? 7 a 


20 MLE 0.0476 0.1718 -0.0669 
BCE -0.0016 -0.0081 -0.0061 


0.6984 0.3947 0.0859 
0.5264 0.2783 0.0847 


30 MLE 0.0313 0.1077 -0.0439 
BCE 0.0004 0.0012 -0.0024 


0.5245 0.2750 0.0684 
0.4478 0.2252 0.0678 


40 MLE 0.0215 0.0754 -0.0330 
BCE -0.0001 -0.0003 -0.0015 


0.4222 0.2207 0.0582 
0.3835 0.1954 0.0578 


50 MLE 0.0160 0.0558 -0.0259 
BCE 0.0000 -0.0001 -0.0005 


0.3609 0.1862 0.0516 
0.3380 0.1710 0.0514 



BCE: bias corrected estimate. 
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Fig. 1. Scatter-plot of the data set. 

6 Application 

Obviously, due to the genesis of the Birnbaum-Saunders distribution, the 
fatigue processes are by excellence ideally modeled by this model. We now 
consider an application to a biaxial fatigue data set reported by Ricck and 
Nedelman (1991) on the life of a metal piece in cycles to failure. The response 
N is the number of cycles to failure and the explanatory variable w is the 
work per cycle (mJ/m^). The data of forty six observations were taken from 
Table 1 of Galea et al. (2004). 

Rieck and Nedelman (1991) proposed the following model for the biaxial fa- 
tigue data: 

= /3i + /32logWi + £j, (9) 
where = logA^^ and £j ~ W(a,0,2), for i = 1,...,46. The MLEs (the 
corresponding standard errors in parentheses) are: f3i = 12.2797 (0.3942), /?2 = 
-1.6708(0.1096) and a = 0.4104(0.0428). We take the logarithm of w to 
ensure a linear relationship between the response variable (y) and the covariate 
in (9); see Galea et al. (2004, Figure 1). However, Figure 1 suggests a nonlinear 
relationship between the response variable and the covariate w. 

Here, we proposed the nonlinear regression model 

= /3i + /?2exp(/33M)+£i, i = l,...,46, (10) 

where Si ~ SN'{a, 0,2). The MLEs (the standard errors in parentheses) are: 
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Fig. 2. Scatter-plot and the fitted models. 

Pi = 8.9876(0.7454), P2 = -5.1802(0.5075), = -22.5196 (7.3778) jind 
a = 0.40 (0.0417). The bias corrected estimates are: Pi = 8.7806 (0.7734), p2 = 
-4.9362 (0.5266), P3 = -22.1713 (7.6548) and a = 0.4157 (0.0433). Hence, the 
uncorrected estimates are shghtly different from the bias corrected estimates 
even for large samples (n = 46 observations). 

Figure 2 gives the scatter-plot of the data, the fitted model (10) and the fitted 
straight line, say yi = (3i+f32Wi+ei, where the MLEs are: pi = 7.9864 (0.1622), 
p2 = -0.0406(0.0036) and a = 0.52(0.0542). Figure 2 shows that the non- 
linear model (10) (unlike the linear model) fits satisfactorily to the fatigue 
data. The 46th observation (the one with work per cycle near 100) can be an 
influential data. However, it is not possible to say whether this observation is 
influential or not without using an efficient way to detect influential observa- 
tions in the new class of models. Influence diagnostic analysis for this class of 
models will be developed in future research. 

Following Xie and Wei (2007), we obtain the residuals Si = yi — fii and 
Ri — 2S~^sinh(ej/2). Figure 3 gives the scatter-plot of Ri versus the pre- 
dicted values pi for both fltted models: (i) yi — Pi + P2Wi + ef, and (ii) 
Hi — Pi P2'^^p{P3/wi) + Ei- Figure 3 shows that the distribution of Ri is 
approximately normal for model (ii) but this is not true for model (i). Based 
upon the fact that U ~ SJ\f{a, fi, a) if 2a~^ sinh{(t/ — ii)/a} ~ A/'(0, 1), then 
the residual ii should follow approximately a sinh-normal distribution. 
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Fig. 3. Index plot of Ri versus juj. 

7 Conclusions 

The Birnbaum-Saunders distribution is widely used to model times to fai- 
lure for materials subject to fatigue. The purpose of the paper was two fold. 
First, we propose a new class of Birnbaum-Saunders nonlinear regression mo- 
dels which generalizes the regression model described in Rieck and Nedelman 
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(1991). Second, wc give simple formulae for calculating bias corrected ma- 
ximum likelihood estimates of the parameters of these models. The simulation 
results presented show that the bias correction derived is very effective, even 
when the sample size is large. Indeed, the bias correction mechanism adopted 
yields adjusted maximum likelihood estimates which are nearly unbiased. We 
also present an application to a real fatigue data set that illustrates the use- 
fulness of the proposed model. Future research will be devoted to a study of 
diagnostics and influence analysis in the new class of nonlinear models. 
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